
clear
set more off

**************************************************
*PREPARE DATA
**************************************************

insheet using "output-base/identification_mcaidpref.csv"
xpose, clear
rename v1 x1
rename v2 y1
rename v3 y2
label var x1 "Medicaid preference"
label var y1 "Informal care rate, Q1"
label var y2 "Informal care rate, Q2"
sort x1
tempfile p1
save `p1', replace

insheet using "output-base/identification_altruism.csv", clear
xpose, clear
rename v1 x2
rename v2 y2_altruism
rename v4 y22
rename v6 y23
label var x2 "Altruism"
label var y2_altruism "Median parent assets, age 79"
sort x2
tempfile p2
save `p2', replace

insheet using "output-base/identification_stigma.csv", clear
xpose, clear
sort v1
rename v1 x3
rename v7 y3
rename v6 y3_2
label var x3 "Child guilt"
label var y3 "LTC insurance rate, avg"
label var y3_2 "LTC insurance rate, Q5"
sort x3
tempfile p3
save `p3', replace


insheet using "output-base/identification_formalpref.csv", clear
xpose, clear
sort v1
rename v1 x4
rename v4 y4_3
rename v7 y4_avg
label var x4 "Formal care preference"
sort x4
tempfile p4
save `p4', replace

insheet using "output-base/identification_pareto.csv", clear
xpose, clear
sort v1
rename v1 x5
rename v6 y5
label var x5 "output/Initial pareto weight"
sort x5
tempfile p5
save `p5', replace

insheet using "output-base/identification_altruism_ratiomoments.csv", clear
xpose, clear
sort v1
rename v1 x6
rename v6 y6
label var x6 "Altruism"
label var y6 "Ratio"
sort x6
tempfile p6
save `p6', replace

insheet using "output-base/identification_pareto_wealthmoments.csv", clear
xpose, clear
sort v1
rename v1 x7
rename v2 y7
label var x7 "Initial pareto weight"
label var y7 "Wealth"
sort x7
tempfile p7
save `p7', replace

insheet using "output-base/identification_paretoaltruismjoint.csv", clear
xpose, clear
sort v1
rename v1 x8
rename v2 y8
label var x8 "Initial pareto weight"
sort x8
tempfile p8
save `p8', replace


**************************************************


**************************************************
*MERGE DATA AND CREATE GRAPHS!
**************************************************
use `p1', clear
merge 1:1 _n using `p4', nogen
merge 1:1 _n using `p3', nogen
merge 1:1 _n using `p2', nogen
merge 1:1 _n using `p5', nogen
merge 1:1 _n using `p6', nogen
merge 1:1 _n using `p7', nogen
merge 1:1 _n using `p8', nogen

cap erase p1.gph
cap erase p3.gph
cap erase p4.gph
cap erase p2.gph
cap erase p5.gph
cap erase p6.gph
cap erase p7.gph
cap erase p8.gph

*MEDICAID PREFERENCE
gen data1=.3927813
gen est1=-44.5E-007
replace est1=est1/1E-7
replace x1=x1/1E-7
twoway line data1 x1, connect(direct) lpattern(dash) lwidth(thick) lcolor(lavender) ///
	|| line y1 x1, connect(direct) lwidth(thick) lcolor(black) ///
	|| line y1 est1, connect(direct) lwidth(thick) lcolor(lavender) ///
	   ytitle("Informal care rate, Q1") xtitle(Medicaid preference (*10{superscript:-7})) ///
	   xlabel(-100(25)25) saving(p1.gph) ///
	   graphregion(color(white)) legend(label(1 "Data moment") label(2 "Simulated moment") label(3 "Estimate") rows(1))
	   
*ALTRUISM  
gen data2=78919.58/1000
replace y2_altruism=y2_altruism/1000
local est2=7.86604573807630075E-002
twoway line y2_altruism x2, lwidth(thick) lcolor(black) 	|| line data2 x2,  lpattern(dash) lwidth(thick) lcolor(lavender) ///
	   xline(`est2', lwidth(thick) lcolor(lavender)) ///
	   ytitle(Median parent wealth at age 79) xtitle(Parent altruism) ///
	   saving(p2.gph) graphregion(color(white)) legend(off)

*STIGMA	 
gen data3=.0769298
gen est3=3E-7/1E-7
replace x3=-x3/1E-7
twoway line y3 x3, lwidth(thick) lcolor(black) ///
	|| line data3 x3, connect(direct) lpattern(dash) lwidth(thick) lcolor(lavender) ///
	|| line y3 est3, connect(direct) lwidth(thick) lcolor(lavender)  ///
	   xlabel(-5(4)9) xtitle(Child guilt (*10{superscript:-7}))  ytitle("Insurance rate, average")  ///
	   saving(p3.gph) graphregion(color(white)) legend(label(1 "Simulated moment") label(2 "Estimate") rows(1))

*FORMAL CARE PREF
gen data4_avg=.4982489
gen data4_3=.6164384
replace x4=x4/1E-7
local est4=-2E-7/1E-7
twoway line data4_3 x4, lpattern(dash) lwidth(thick) lcolor(lavender) ///
	|| line y4_3 x4, lwidth(thick) lcolor(black) ///
	|| line data4_avg x4, lpattern(dash) lwidth(thick) lcolor(lavender*.5) ///
	|| line y4_avg x4,  lwidth(thick) lcolor(gs7) ///
	xline(`est4', lwidth(thick) lcolor(lavender)) ///
   ytitle("Informal care rate (median, avg)") xtitle(Formal care preference (*10{superscript:-7})) ///
   graphregion(color(white)) saving(p4.gph)
   
*INITIAL PARETO WEIGHT
gen data5=103553.8/60020.27
local est5=.6
recode x5 (5=.4) (6=.5) (7=.6) (8=.7) (9=.8)
twoway line data5 x5, lpattern(dash) lwidth(thick) lcolor(lavender) ///
	|| line y5 x5, lwidth(thick) lcolor(black) ///
	xline(`est5', lwidth(thick) lcolor(lavender)) ///
   ytitle("Ratio of parent to child assets, age 71") xtitle("Initial Pareto weight") ///
   graphregion(color(white)) saving(p5.gph)   
   
cap graph drop toresize
grc1leg p1.gph p2.gph p3.gph p4.gph p5.gph, graphregion(color(white)) name(toresize) iscale(.75) row(2) legendfrom(p1.gph)
graph display toresize //, ysize(5) xsize(5)
graph export "identification_graphs.pdf", as(pdf) replace 

*ALTRUISM ratio	   
gen data6=103553.8/60020.27
gen est6=7.86604573807630075E-002
twoway line y6 x6, lwidth(thick) lcolor(black) ||	///
	   line data6 x6,  lpattern(dash) lwidth(thick) lcolor(lavender) || ///
	   line y6 est6, connect(direct) lwidth(thick) lcolor(lavender) ///
	   ytitle("Ratio of parent to child assets, age 71") xtitle(Parent altruism) ///
	   legend(label(1 "Data moment") label(2 "Simulated moment") label(3 "Estimate") rows(1)) ///
	   graphregion(color(white))
graph export "identification_graphs_altruismratio.pdf", as(pdf) replace
  
*Pareto wealth  
gen data7=78919.58/1000
replace y7=y7/1000
gen est7=.6
recode x7 (5=.4) (6=.5) (7=.6) (8=.7) (9=.8)
twoway line y7 x7, lwidth(thick) lcolor(black) 	|| ///
	   line data7 x7,  lpattern(dash) lwidth(thick) lcolor(lavender) || ///
	   line y7 est7, connect(direct) lwidth(thick) lcolor(lavender) ///
	   ytitle(Median parent wealth at age 79) xtitle("Initial Pareto weight") ///
	   legend(label(1 "Data moment") label(2 "Simulated moment") label(3 "Estimate") rows(1)) ///
	   graphregion(color(white)) 
graph export "identification_graphs_paretowealth.pdf", as(pdf) replace
 
cap graph drop toresize
grc1leg p6.gph p7.gph, graphregion(color(white)) name(toresize) iscale(.75) row(2) legendfrom(p6.gph)
graph display toresize //, ysize(5) xsize(5)
graph export "identification_graphs_paretoaltruism.pdf", as(pdf) replace  
 
*INITIAL PARETO WEIGHT, TILDE (JOINT WITH ALTRUISM)
gen data8=103553.8/60020.27
gen est8=.6
recode x8 (7=.56) (8=.58) (9=.6) (10=.62) (11=.64)
twoway line data8 x8, lpattern(dash) lwidth(thick) lcolor(lavender) ///
	|| line y8 x8, lwidth(thick) lcolor(black) ///
	|| line y8 est8, connect(direct) lwidth(thick) lcolor(lavender) ///
   ytitle("Ratio of parent to child assets, age 71") xtitle("Initial Pareto weight (altruism adjusted)") ///
   legend(label(1 "Data moment") label(2 "Simulated moment") label(3 "Estimate") rows(1)) /// 
   graphregion(color(white))
graph export "identification_graphs_paretoaltruismjoint.pdf", as(pdf) replace

**************************************************

*SENSITIVITY ANALYSIS (ANDREWS GENTZKOW SHAPIRO)
*note: have to print out weight matrix and sensitivity from fortran (just variances), save as csv
insheet using "varmoments.csv", clear nonames
rename v1 varmoments
destring varmoments, replace force
replace varmoments=.0000000139 if _n==1 //transfer issue
tempfile vars
save `vars', replace

insheet using "sensitivity_AGS.csv", clear nonames double
xpose, clear
merge 1:1 _n using `vars', nogen
gen double sdmoments=sqrt(varmoments)

gen 	momenttype=1 if inlist(_n,1,2,3,4,5,10,11)	//kid assets young (<=75)
replace momenttype=2 if inlist(_n,16,17,23,24,25,26,27,29,30,31,32,33,34,35,36,37,38,39) //kid assets (85+)
replace momenttype=3 if inlist(_n,39+1,39+2,39+3,39+4,39+5,39+10,39+11) 	//parent assets young (<=75)
replace momenttype=4 if inlist(_n,39+16,39+17,39+23,39+24,39+25,39+26,39+27, ///
								   39+29,39+30,39+31,39+32,39+33,39+34,39+35,39+36,39+37,39+38,39+39) 	//parent assets old (85+)
replace momenttype=5 if inrange(_n,79,80)	//insurance
replace momenttype=6 if inrange(_n,82,83)	//insurance
replace momenttype=7 if inrange(_n,84,85)	//informal care
replace momenttype=8 if inrange(_n,87,88)	//informal care
replace momenttype=9 if _n==89				//medicaid
replace momenttype=10 if inrange(_n,90,94)	//labor supply

local pctchange1=-2E-7 					   //formal care preference                     
local pctchange2=3E-7   				   //stigma
local pctchange3=7.86604573807630075E-002  //altruism
local pctchange4= 2.26  				   //alpha
local pctchange5=0.54142202266996531       //L shifter
local pctchange6=1.8888369315929308        //gammab
local pctchange7=2.1143038768435756         //gamma
local pctchange8=23381				        //cminp
local pctchange9=-4.45076144931274663E-006 //mcaidpref
local pctchange10=.6 						//pareto 

forval p=1/10 {
	gen double v`p'_pctchange=-((v`p')/`pctchange`p'')/sdmoments
}
forval p=1/10 {
	bys momenttype: egen double momtypeavg`p'=mean(v`p'_pctchange)
}
keep momtype*
duplicates drop

mkmat momtypeavg* if inrange(_n,1,10), mat(sensitivity)
heatplot sensitivity, color(hcl, diverging) cuts(-.002(.0005).002) ///
	xlabel(1 "z{sub:F}" 2 "g" 3 "{&eta}" 4 "{&gamma}{sub:ell}" 5 "{&alpha}" 6 "{&gamma}{sub:b}" ///
		   7 "{&gamma}{sub:c}" 8 "floor" 9 "z{sub:M}" 10 "{&theta}{sub:0}", angle(45)) ///
	ylabel(1 "Child assets, young" 2 "Child assets, old" 3 "Parent assets, young" ///
		   4 "Parent assets, old" 5 "Insurance, poor" 6 "Insurance, rich" ///
		   7 "Informal care, poor" 8 "Informal care, rich" 9 "Medicaid" 10 "Full-time") ///
		   ytitle(Moment) xtitle(Parameter) graphregion(color(white))
graph export "identification_graphs_AGSsensitivity.pdf", as(pdf) replace
